library(readr)
package 㤼㸱readr㤼㸲 was built under R version 4.0.5
library(ggplot2)
library(dplyr)
library(rms)
package 㤼㸱rms㤼㸲 was built under R version 4.0.5package 㤼㸱survival㤼㸲 was built under R version 4.0.5package 㤼㸱SparseM㤼㸲 was built under R version 4.0.4
# library(synthpop)
library(tidyr)
package 㤼㸱plotly㤼㸲 was built under R version 4.0.5
library(sjPlot)
package 㤼㸱sjPlot㤼㸲 was built under R version 4.0.5
library(performance)
# install the 'mi' package if you are planning on redoing imputation on the original dataset. For now the imputed dataset has already been saved
# library(mi)
Reading in original dataset, performing multiple imputation and writing out imputed dataset. This chunk commented out because we’ve already done this and saved the resulting file as our starting point.
Read in dataset from Brinati paper. Dataset already been imputed based on code above.
brinati_covid_study_v2.imputed <- read_csv("data/brinati-covid_study_v2_imputed.csv",
col_types = cols(GENDER = col_factor(levels = c("M",
"F")), SWAB = col_factor(levels = c("0",
"1"))))
glm.orig.fit<-glm(SWAB ~ ., brinati_covid_study_v2.imputed, family='binomial')
glm.fit: fitted probabilities numerically 0 or 1 occurred
glm.orig.fit
Call: glm(formula = SWAB ~ ., family = "binomial", data = brinati_covid_study_v2.imputed)
Coefficients:
(Intercept) GENDERF AGE WBC Platelets Neutrophils Lymphocytes
-2.185e+00 -8.505e-01 1.746e-02 -1.146e-01 2.011e-03 -1.289e-01 9.374e-03
Monocytes Eosinophils Basophils CRP AST ALT ALP
8.215e-01 -1.095e-08 -4.520e+00 5.427e-04 -7.082e-03 1.026e-02 -1.195e-02
1.001e-02
Degrees of Freedom: 278 Total (i.e. Null); 263 Residual
Null Deviance: 366.4
Residual Deviance: 239.4 AIC: 271.4
write.csv(summary(glm.orig.fit)$coef, file='output/OR-original-model.csv', row.names = FALSE)
predicted = plogis(predict(glm.orig.fit))
C (ROC) R2 D D:Chi-sq D:p U
7.322477e-01 8.661239e-01 5.000567e-01 4.514185e-01 1.269458e+02 NA -7.168459e-03
U:Chi-sq U:p Q Brier Intercept Slope Emax
7.673862e-13 1.000000e+00 4.585870e-01 1.361516e-01 -2.780966e-08 1.000000e+00 6.627578e-02
E90 Eavg S:z S:p
5.261796e-02 2.628540e-02 -2.750989e-01 7.832403e-01
p<-plot_model(glm.orig.fit)
Argument 'df_method' is deprecated. Please use 'ci_method' instead.
p<-ApplyFigureThemeLargeFontOnly(p + ylim(.1, 2.5) )
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing
scale.
p
ggsave(filename = 'figs/OR-original-model.png', plot=p, width=6, height=6, units='in', dpi=300)

Draw ROC curve
cms.cutoffs <- lapply(seq(0.01,0.99,0.01), function(cutoff) {
ret=confusionMatrix(predicted, observed, cutoff=cutoff)
ret$cutoff = cutoff
ret
})
ff<-data.frame(t(sapply(cms.cutoffs, function(cf) c("Cutoff"=cf$cutoff, "Sensitivty"=cf$sens, "Specificity"=cf$spec))))
cutoffs_to_plot <- c(0.09, 0.3, 0.6, 0.9)
geom_point(size=3, color='darkblue') + xlim(0,1) + ylim(0,1)+
xlab('False Positive Rate\n(1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')
)
ggsave(filename = 'figs/example-sens-spec-on-roc.png', plot=p, width=6, height=6, units='in', dpi=300)
for (cutoff in lapply(1:length(cutoffs_to_plot), function(i) cutoffs_to_plot[1:i])) {
p <- ApplyFigureTheme(ff %>% mutate(fpr=1-Specificity) %>% filter(Cutoff %in% cutoff) %>%
ggplot(., aes(x=fpr, y=Sensitivity)) +
geom_point(size=3, color='darkblue') + xlim(0,1) + ylim(0,1)+
xlab('False Positive Rate\n(1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')
)
ggsave(filename = paste0('figs/example-sens-spec-on-roc-', cutoff[length(cutoff)], '.png'),
plot=p, width=6, height=6, units='in', dpi=300)
}
p <- ff %>% mutate(fpr=1-Specificity) %>%
ggplot(., aes(x=fpr, y=Sensitivity)) + geom_point(aes(frame=Cutoff)) +
xlab('False Positive Rate (1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')

Draw calibration
cal_breaks = seq(0, 1, .1)
cal_deciles <- data.frame(Predicted=predicted, Observed = as.integer(as.character(observed))) %>%
mutate(Bins = cut(Predicted, breaks=cal_breaks, include.lowest = TRUE))
cal_deciles_summary <- cal_deciles %>% group_by(Bins) %>% summarise(n=n(), Average_Observed = mean(Observed),
Average_Predicted=mean(Predicted))
for (cutoff in lapply(1:nrow(cal_deciles_summary), function(i) cal_deciles_summary$Bins[1:i])) {
SaveStdSquareFigure(p, filename = paste0('figs/example-cal-curve-points-', cutoff[length(cutoff)], '.png'))
p<-ApplyFigureThemeCalCurvePoints(ggplot(cal_deciles_summary , aes(x=Average_Predicted, y=Average_Observed)) )
p<-ApplyFigureThemeLargeFontOnly(p)
SaveStdSquareFigure(p, 'figs/cal-curve-deciles-all.png')
p<-p + geom_smooth(method='lm', se=FALSE)
SaveStdSquareFigure(p, 'figs/cal-curve-deciles-all-with-bestfit.png')
p

NA
Calculate hosmer lemeshow statistics
hl.org <- performance_hosmer(glm.orig.fit)
hl.org
# Hosmer-Lemeshow Goodness-of-Fit Test
Chi-squared: 5.333
df: 8
p-value: 0.722
Summary: model seems to fit well.
Create smoothed cal curve from Van hoorde et al
p<-CreateSmoothedCalCurvePlot(predicted, observed)
p<-ApplyFigureThemeLargeFontOnly(p)
SaveStdSquareFigure(p, 'figs/smoothed-cal-curve-orig-data.png')
p


Boot strap performance in the original dataset
Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls

Assess performance in external validation
- Look at performance as a function of prevalence
prevalences = seq(0.01, 0.99, 0.01)
brinati.syn.factory = generateSyntheticDataFactory(baseDF = brinati_covid_study_v2.imputed, method = 'boot')
dfs_prevs <- lapply(prevalences, function (prev) brinati.syn.factory(prev = prev))
dfs.prevs.perf <- lapply(dfs_prevs, function(df) {
assessPerf(predicted = plogis(predict(glm.orig.fit, newdata=as.data.frame(df))), observed = df$SWAB)
})
Calculate the different validation measures
dfs.prevs.perf.df <- data.frame(matrix(unlist(dfs.prevs.perf), nrow = length(prevalences), byrow = TRUE))
colnames(dfs.prevs.perf.df) <- names(dfs.prevs.perf[[1]])
dfs.prevs.perf.df$Prevalence = prevalences
dfs.prevs.perf.df <- gather(dfs.prevs.perf.df, Measure, Value, Dxy:`S:p`, factor_key = TRUE)
Plot the distribution of validation measures
measures = c('C (ROC)', 'Brier', 'Intercept', 'Slope')
p<-ggplot(dfs.prevs.perf.df %>% filter(Measure %in% measures), aes(x=Measure, y=Value, color=Measure)) + geom_boxplot() +
geom_point(data=data.frame(Measure = measures, Value = c(1, 0, 0, 1)), size=3, color='blue') + theme_bw()
p<-ApplyFigureThemeLargeFontOnly(p)
SaveMedRectFigure(p, 'figs/change-in-metrics-due-to-prevalence-boxplot.png')
p
p<-ggplot(dfs.prevs.perf.df %>% filter(Measure %in% measures), aes(x=Prevalence, y=Value, color=Measure, group=Measure)) + geom_line(size=1.3) + theme_bw()
p<-ApplyFigureThemeLargeFontOnly(p)
SaveMedRectFigure(p, 'figs/change-in-metrics-due-to-prevalence-linearplot.png')
p
Calculate Decision Curve
cms.cutoffs <- lapply(seq(0.01,0.99,0.01), function(cutoff) {
ret=confusionMatrix(predicted, observed, cutoff=cutoff)
ret$cutoff = cutoff
ret
})
calculate_net_benefit_at_threshold_for_sensitivity_and_specificity <- function(sens, spec, threshold) {
sens - (1-spec)*(threshold/(1-threshold))
}
netbenefit_across_thresholds = sapply(cms.cutoffs, function(cutoff_values) calculate_net_benefit_at_threshold_for_sensitivity_and_specificity(cutoff_values$sens, cutoff_values$spec, cutoff_values$cutoff))
netbenefit_vs_thresholds = data.frame(thresholds = sapply(cms.cutoffs, function(cutoff_values) cutoff_values$cutoff), netbenefit = netbenefit_across_thresholds)
ggplot(netbenefit_vs_thresholds, aes(x=thresholds, y=netbenefit)) + geom_line() + ylab("Net Benefit") + xlab("Threshold to Act") + theme_bw()
LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KYGBge3IgbWVzc2FnZT1GQUxTRX0NCmxpYnJhcnkocmVhZHIpDQpsaWJyYXJ5KGdncGxvdDIpDQpsaWJyYXJ5KGRwbHlyKQ0KbGlicmFyeShybXMpDQojIGxpYnJhcnkoc3ludGhwb3ApDQpsaWJyYXJ5KHRpZHlyKQ0KbGlicmFyeShwbG90bHkpDQpsaWJyYXJ5KHNqUGxvdCkNCmxpYnJhcnkocGVyZm9ybWFuY2UpDQojIGluc3RhbGwgdGhlICdtaScgcGFja2FnZSBpZiB5b3UgYXJlIHBsYW5uaW5nIG9uIHJlZG9pbmcgaW1wdXRhdGlvbiBvbiB0aGUgb3JpZ2luYWwgZGF0YXNldC4gRm9yIG5vdyB0aGUgaW1wdXRlZCBkYXRhc2V0IGhhcyBhbHJlYWR5IGJlZW4gc2F2ZWQNCiMgbGlicmFyeShtaSkNCnNvdXJjZSgnLi9oZWxwZXIgZnVuY3Rpb25zLlInKQ0KYGBgDQoNClJlYWRpbmcgaW4gb3JpZ2luYWwgZGF0YXNldCwgcGVyZm9ybWluZyBtdWx0aXBsZSBpbXB1dGF0aW9uIGFuZCB3cml0aW5nIG91dCBpbXB1dGVkIGRhdGFzZXQuIFRoaXMgY2h1bmsgY29tbWVudGVkIG91dCBiZWNhdXNlIHdlJ3ZlIGFscmVhZHkgZG9uZSB0aGlzIGFuZCBzYXZlZCB0aGUgcmVzdWx0aW5nIGZpbGUgYXMgb3VyIHN0YXJ0aW5nIHBvaW50Lg0KYGBge3IgZWNobz1GQUxTRSwgd2FybmluZz1GQUxTRX0NCiMgDQojIGJyaW5hdGlfY292aWRfc3R1ZHlfdjIgPC0gcmVhZF9jc3YoImRhdGEvYnJpbmF0aS1jb3ZpZF9zdHVkeV92Mi5jc3YiLA0KIyBjb2xfdHlwZXMgPSBjb2xzKEdFTkRFUiA9IGNvbF9mYWN0b3IobGV2ZWxzID0gYygiTSIsDQojICJGIikpLCBTV0FCID0gY29sX2ZhY3RvcihsZXZlbHMgPSBjKCIwIiwNCiMgIjEiKSkpKQ0KIyBicmluYXRpX2NvdmlkX3N0dWR5X3YyIDwtIGJyaW5hdGlfY292aWRfc3R1ZHlfdjIgJT4lIG11dGF0ZV9pZihpcy5udW1lcmljLCAuZnVucyA9IGZ1bmN0aW9uKHgpIHt4KzAuMDAxfSkNCiMgYnJhaW5hdGkubWkgPC0gbWlzc2luZ19kYXRhLmZyYW1lKGFzLmRhdGEuZnJhbWUoYnJpbmF0aV9jb3ZpZF9zdHVkeV92MiksIGZhdm9yX3Bvc2l0aXZlID0gVFJVRSkNCiMgDQojIA0KIyBicmFpbmF0aS5taTwtY2hhbmdlKGJyYWluYXRpLm1pLCB5PWMoJ0x5bXBob2N5dGVzJywgJ0Jhc29waGlscycsICdNb25vY3l0ZXMnLCdFb3Npbm9waGlscycsICdCYXNvcGhpbHMnKSwgd2hhdD0ndHlwZScsIHRvPSdwb3NpdGl2ZS1jb250aW51b3VzJykNCiMgc2hvdyhicmFpbmF0aS5taSkNCiMgaW1hZ2UoYnJhaW5hdGkubWkpDQojIA0KIyBvcHRpb25zKG1jLmNvcmVzID0gNCkNCiMgaW1wdXRhdGlvbnMgPC0gbWkoYnJhaW5hdGkubWksIG4uaXRlciA9IDkwLCBuLmNoYWlucyA9IDQsIG1heC5taW51dGVzID0gMjApIA0KIyBzaG93KGltcHV0YXRpb25zKQ0KIyByb3VuZChtaXBwbHkoaW1wdXRhdGlvbnMsIG1lYW4sIHRvLm1hdHJpeCA9IFRSVUUpLCAzKQ0KIyBSaGF0cyhpbXB1dGF0aW9ucykNCiMgZGZzIDwtIGNvbXBsZXRlKGltcHV0YXRpb25zLCBtPTIpDQojIGJyaW5hdGlfY292aWRfdjIuaW1wdXRlZCA8LSBkZnNbWzFdXVssMToxNl0NCiMgd3JpdGUuY3N2KGJyaW5hdGlfY292aWRfdjIuaW1wdXRlZCwgZmlsZT0nLi9kYXRhL2JyaW5hdGktY292aWRfc3R1ZHlfdjJfaW1wdXRlZC5jc3YnLCByb3cubmFtZXMgPSBGQUxTRSkNCmBgYA0KDQpSZWFkIGluIGRhdGFzZXQgZnJvbSBCcmluYXRpIHBhcGVyLiBEYXRhc2V0IGFscmVhZHkgYmVlbiBpbXB1dGVkIGJhc2VkIG9uIGNvZGUgYWJvdmUuIA0KYGBge3J9DQoNCmJyaW5hdGlfY292aWRfc3R1ZHlfdjIuaW1wdXRlZCA8LSByZWFkX2NzdigiZGF0YS9icmluYXRpLWNvdmlkX3N0dWR5X3YyX2ltcHV0ZWQuY3N2IiwNCiAgY29sX3R5cGVzID0gY29scyhHRU5ERVIgPSBjb2xfZmFjdG9yKGxldmVscyA9IGMoIk0iLA0KICAiRiIpKSwgU1dBQiA9IGNvbF9mYWN0b3IobGV2ZWxzID0gYygiMCIsDQogICIxIikpKSkNCmBgYA0KDQoNCmBgYHtyfQ0KDQpnbG0ub3JpZy5maXQ8LWdsbShTV0FCIH4gLiwgYnJpbmF0aV9jb3ZpZF9zdHVkeV92Mi5pbXB1dGVkLCBmYW1pbHk9J2Jpbm9taWFsJykNCmdsbS5vcmlnLmZpdA0Kd3JpdGUuY3N2KHN1bW1hcnkoZ2xtLm9yaWcuZml0KSRjb2VmLCBmaWxlPSdvdXRwdXQvT1Itb3JpZ2luYWwtbW9kZWwuY3N2Jywgcm93Lm5hbWVzID0gRkFMU0UpDQpwcmVkaWN0ZWQgPSBwbG9naXMocHJlZGljdChnbG0ub3JpZy5maXQpKQ0Kb2JzZXJ2ZWQgPSBicmluYXRpX2NvdmlkX3N0dWR5X3YyLmltcHV0ZWQkU1dBQg0KYXNzZXNzUGVyZihwcmVkaWN0ZWQsIG9ic2VydmVkKQ0KcDwtcGxvdF9tb2RlbChnbG0ub3JpZy5maXQpDQpwPC1BcHBseUZpZ3VyZVRoZW1lTGFyZ2VGb250T25seShwICsgeWxpbSguMSwgMi41KSApIA0KICANCnANCmdnc2F2ZShmaWxlbmFtZSA9ICdmaWdzL09SLW9yaWdpbmFsLW1vZGVsLnBuZycsIHBsb3Q9cCwgd2lkdGg9NiwgaGVpZ2h0PTYsIHVuaXRzPSdpbicsIGRwaT0zMDApDQpgYGANCg0KIyBEcmF3IFJPQyBjdXJ2ZQ0KYGBge3J9DQpjbXMuY3V0b2ZmcyA8LSBsYXBwbHkoc2VxKDAuMDEsMC45OSwwLjAxKSwgZnVuY3Rpb24oY3V0b2ZmKSB7DQogIHJldD1jb25mdXNpb25NYXRyaXgocHJlZGljdGVkLCBvYnNlcnZlZCwgY3V0b2ZmPWN1dG9mZikNCiAgcmV0JGN1dG9mZiA9IGN1dG9mZg0KICByZXQNCiAgfSkNCmZmPC1kYXRhLmZyYW1lKHQoc2FwcGx5KGNtcy5jdXRvZmZzLCBmdW5jdGlvbihjZikgYygiQ3V0b2ZmIj1jZiRjdXRvZmYsICJTZW5zaXRpdnR5Ij1jZiRzZW5zLCAiU3BlY2lmaWNpdHkiPWNmJHNwZWMpKSkpDQpjb2xuYW1lcyhmZikgPC0gYygnQ3V0b2ZmJywgJ1NlbnNpdGl2aXR5JywnU3BlY2lmaWNpdHknKQ0KDQoNCmN1dG9mZnNfdG9fcGxvdCA8LSBjKDAuMDksIDAuMywgMC42LCAwLjkpDQpwIDwtIEFwcGx5RmlndXJlVGhlbWUoZmYgJT4lIG11dGF0ZShmcHI9MS1TcGVjaWZpY2l0eSkgJT4lIGZpbHRlcihDdXRvZmYgJWluJSBjdXRvZmZzX3RvX3Bsb3QpICU+JQ0KICBnZ3Bsb3QoLiwgYWVzKHg9ZnByLCB5PVNlbnNpdGl2aXR5KSkgKyANCiAgICBnZW9tX3BvaW50KHNpemU9MywgY29sb3I9J2RhcmtibHVlJykgKyB4bGltKDAsMSkgKyB5bGltKDAsMSkrDQogICAgeGxhYignRmFsc2UgUG9zaXRpdmUgUmF0ZVxuKDEtU3BlY2lmaWNpdHkpJykgKyB5bGFiKCdUcnVlIFBvc2l0aXZlIFJhdGVcbihTZW5zaXRpdml0eSknKQ0KICApDQpnZ3NhdmUoZmlsZW5hbWUgPSAnZmlncy9leGFtcGxlLXNlbnMtc3BlYy1vbi1yb2MucG5nJywgcGxvdD1wLCB3aWR0aD02LCBoZWlnaHQ9NiwgdW5pdHM9J2luJywgZHBpPTMwMCkNCg0KZm9yIChjdXRvZmYgaW4gbGFwcGx5KDE6bGVuZ3RoKGN1dG9mZnNfdG9fcGxvdCksIGZ1bmN0aW9uKGkpIGN1dG9mZnNfdG9fcGxvdFsxOmldKSkgew0KICBwIDwtIEFwcGx5RmlndXJlVGhlbWUoZmYgJT4lIG11dGF0ZShmcHI9MS1TcGVjaWZpY2l0eSkgJT4lIGZpbHRlcihDdXRvZmYgJWluJSBjdXRvZmYpICU+JQ0KICAgIGdncGxvdCguLCBhZXMoeD1mcHIsIHk9U2Vuc2l0aXZpdHkpKSArIA0KICAgICAgZ2VvbV9wb2ludChzaXplPTMsIGNvbG9yPSdkYXJrYmx1ZScpICsgeGxpbSgwLDEpICsgeWxpbSgwLDEpKw0KICAgICAgeGxhYignRmFsc2UgUG9zaXRpdmUgUmF0ZVxuKDEtU3BlY2lmaWNpdHkpJykgKyB5bGFiKCdUcnVlIFBvc2l0aXZlIFJhdGVcbihTZW5zaXRpdml0eSknKQ0KICApDQogIGdnc2F2ZShmaWxlbmFtZSA9IHBhc3RlMCgnZmlncy9leGFtcGxlLXNlbnMtc3BlYy1vbi1yb2MtJywgY3V0b2ZmW2xlbmd0aChjdXRvZmYpXSwgICcucG5nJyksIA0KICAgICAgICAgcGxvdD1wLCB3aWR0aD02LCBoZWlnaHQ9NiwgdW5pdHM9J2luJywgZHBpPTMwMCkNCn0NCnAgPC0gQXBwbHlGaWd1cmVUaGVtZSgNCiAgZ2dwbG90KGZmLCBhZXMoeD0xLVNwZWNpZmljaXR5LCB5PVNlbnNpdGl2aXR5KSkgKyBnZW9tX2xpbmUoKSsgDQogICAgZ2VvbV9yaWJib24oYWVzKHltaW4gPSAwLCB5bWF4PVNlbnNpdGl2aXR5KSwgY29sb3I9TkEsIGZpbGw9J2JsdWUnLGFscGhhPTAuMykgKyANCiAgICB4bGltKDAsMSkgKyB5bGltKDAsMSkrDQogICAgICB4bGFiKCdGYWxzZSBQb3NpdGl2ZSBSYXRlXG4oMS1TcGVjaWZpY2l0eSknKSArIHlsYWIoJ1RydWUgUG9zaXRpdmUgUmF0ZVxuKFNlbnNpdGl2aXR5KScpIA0KKQ0KZ2dzYXZlKGZpbGVuYW1lID0gJ2ZpZ3MvZXhhbXBsZS1yb2MucG5nJywgDQogICAgICAgICBwbG90PXAsIHdpZHRoPTYsIGhlaWdodD02LCB1bml0cz0naW4nLCBkcGk9MzAwKQ0KcCA8LSBmZiAlPiUgbXV0YXRlKGZwcj0xLVNwZWNpZmljaXR5KSAlPiUNCiAgZ2dwbG90KC4sIGFlcyh4PWZwciwgeT1TZW5zaXRpdml0eSkpICsgZ2VvbV9wb2ludChhZXMoZnJhbWU9Q3V0b2ZmKSkgKw0KICAgIHhsYWIoJ0ZhbHNlIFBvc2l0aXZlIFJhdGUgKDEtU3BlY2lmaWNpdHkpJykgKyB5bGFiKCdUcnVlIFBvc2l0aXZlIFJhdGVcbihTZW5zaXRpdml0eSknKQ0KICANCiAgDQpnZ3Bsb3RseShwKQ0KYGBgDQoNCiMgRHJhdyBjYWxpYnJhdGlvbg0KYGBge3J9DQpjYWxfYnJlYWtzID0gc2VxKDAsIDEsIC4xKQ0KY2FsX2RlY2lsZXMgPC0gZGF0YS5mcmFtZShQcmVkaWN0ZWQ9cHJlZGljdGVkLCBPYnNlcnZlZCA9IGFzLmludGVnZXIoYXMuY2hhcmFjdGVyKG9ic2VydmVkKSkpICU+JQ0KICBtdXRhdGUoQmlucyA9IGN1dChQcmVkaWN0ZWQsIGJyZWFrcz1jYWxfYnJlYWtzLCBpbmNsdWRlLmxvd2VzdCA9IFRSVUUpKQ0KY2FsX2RlY2lsZXNfc3VtbWFyeSA8LSBjYWxfZGVjaWxlcyAlPiUgZ3JvdXBfYnkoQmlucykgJT4lIHN1bW1hcmlzZShuPW4oKSwgQXZlcmFnZV9PYnNlcnZlZCA9IG1lYW4oT2JzZXJ2ZWQpLA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBBdmVyYWdlX1ByZWRpY3RlZD1tZWFuKFByZWRpY3RlZCkpIA0KDQoNCmZvciAoY3V0b2ZmIGluIGxhcHBseSgxOm5yb3coY2FsX2RlY2lsZXNfc3VtbWFyeSksIGZ1bmN0aW9uKGkpIGNhbF9kZWNpbGVzX3N1bW1hcnkkQmluc1sxOmldKSkgew0KICBwPC1BcHBseUZpZ3VyZVRoZW1lQ2FsQ3VydmVQb2ludHMoDQogICAgZ2dwbG90KGNhbF9kZWNpbGVzX3N1bW1hcnkgJT4lIGZpbHRlcihCaW5zICVpbiUgY3V0b2ZmKSwgYWVzKHg9QXZlcmFnZV9QcmVkaWN0ZWQsIHk9QXZlcmFnZV9PYnNlcnZlZCkpICkNCiAgcDwtQXBwbHlGaWd1cmVUaGVtZUxhcmdlRm9udE9ubHkocCkNCiAgU2F2ZVN0ZFNxdWFyZUZpZ3VyZShwLCBmaWxlbmFtZSA9IHBhc3RlMCgnZmlncy9leGFtcGxlLWNhbC1jdXJ2ZS1wb2ludHMtJywgY3V0b2ZmW2xlbmd0aChjdXRvZmYpXSwgICcucG5nJykpDQp9DQoNCnA8LUFwcGx5RmlndXJlVGhlbWVDYWxDdXJ2ZVBvaW50cyhnZ3Bsb3QoY2FsX2RlY2lsZXNfc3VtbWFyeSAsIGFlcyh4PUF2ZXJhZ2VfUHJlZGljdGVkLCB5PUF2ZXJhZ2VfT2JzZXJ2ZWQpKSApDQpwPC1BcHBseUZpZ3VyZVRoZW1lTGFyZ2VGb250T25seShwKQ0KDQpTYXZlU3RkU3F1YXJlRmlndXJlKHAsICdmaWdzL2NhbC1jdXJ2ZS1kZWNpbGVzLWFsbC5wbmcnKQ0KDQpwPC1wICsgZ2VvbV9zbW9vdGgobWV0aG9kPSdsbScsIHNlPUZBTFNFKQ0KU2F2ZVN0ZFNxdWFyZUZpZ3VyZShwLCAnZmlncy9jYWwtY3VydmUtZGVjaWxlcy1hbGwtd2l0aC1iZXN0Zml0LnBuZycpDQpwDQogICAgICAgICAgICANCmBgYA0KDQpDYWxjdWxhdGUgaG9zbWVyIGxlbWVzaG93IHN0YXRpc3RpY3MNCmBgYHtyfQ0KDQpobC5vcmcgPC0gcGVyZm9ybWFuY2VfaG9zbWVyKGdsbS5vcmlnLmZpdCkNCmhsLm9yZw0KYGBgDQoNCkNyZWF0ZSBzbW9vdGhlZCBjYWwgY3VydmUgZnJvbSBWYW4gaG9vcmRlIGV0IGFsDQpgYGB7cn0NCnA8LUNyZWF0ZVNtb290aGVkQ2FsQ3VydmVQbG90KHByZWRpY3RlZCwgb2JzZXJ2ZWQpDQpwPC1BcHBseUZpZ3VyZVRoZW1lTGFyZ2VGb250T25seShwKQ0KU2F2ZVN0ZFNxdWFyZUZpZ3VyZShwLCAnZmlncy9zbW9vdGhlZC1jYWwtY3VydmUtb3JpZy1kYXRhLnBuZycpDQpwDQpgYGANCg0KQm9vdCBzdHJhcCBwZXJmb3JtYW5jZSBpbiB0aGUgb3JpZ2luYWwgZGF0YXNldA0KYGBge3Igd2FybmluZz1GQUxTRSwgZWNobz1GQUxTRX0NCmRmIDwtIGFzLmRhdGEuZnJhbWUoYnJpbmF0aV9jb3ZpZF9zdHVkeV92Mi5pbXB1dGVkKQ0KDQojIHdlIGFyZSBvbmx5IGJvb3RzdHJhcHBpbmcgMjUgdGltZXMgdG8gc2F2ZSB0aW1lOyBpbiByZWFsaXR5IHlvdSB3b3VsZCBpZGVhbGx5IGJvb3RzdHJhcCBhIGZldyBodW5kcmVkIHRpbWVzIGF0IGxlYXN0DQpuQm9vdCA9IDI1DQpzZXQuc2VlZCgxMzQyMzQ5KQ0KdHJhaW5JZHhCb290IDwtIGxhcHBseSgxOm5Cb290LCBmdW5jdGlvbihpKSBzYW1wbGUoMTpucm93KGRmKSwgc2l6ZT1ucm93KGRmKSwgcmVwbGFjZT1UUlVFKSkNCg0KYm9vdFBlcmYub3JpZ2RhdGEgPC0gbGFwcGx5KHRyYWluSWR4Qm9vdCwgDQogICAgICAgICAgICAgICAgICAgZnVuY3Rpb24odHJhaW5JZHgpIA0KICAgICAgICAgICAgICAgICAgICAgYXNzZXNzU2luZ2xlVHJhaW5UZXN0KHRyYWluSWR4LCAoMTpucm93KGRmKSlbLXRyYWluSWR4XSwgZGYsIGdsbS5tb2RlbCA9IFNXQUJ+IC4sIG91dGNvbWUudmFyLm5hbWUgPSAnU1dBQicpDQogICAgICAgICAgICAgICAgICAgKSANCg0KbWVhc3VyZXMgPSBjKCdDIChST0MpJywgJ0JyaWVyJywgJ0ludGVyY2VwdCcsICdTbG9wZScpDQpicmluYXRpLm9yaWcuYm9vdC5wZXJmIDwtIGRhdGEuZnJhbWUobWF0cml4KA0KICB1bmxpc3QobGFwcGx5KGJvb3RQZXJmLm9yaWdkYXRhLCBmdW5jdGlvbihib290KSBib290JHRlc3QucGVyZikpLCANCiAgbnJvdyA9IG5Cb290LCBieXJvdyA9IFRSVUUpKQ0KY29sbmFtZXMoYnJpbmF0aS5vcmlnLmJvb3QucGVyZikgPC0gbmFtZXMoYm9vdFBlcmYub3JpZ2RhdGFbWzFdXSR0ZXN0LnBlcmYpDQoNCmJyaW5hdGkub3JpZy5ib290LnBlcmYgPC0gZ2F0aGVyKGJyaW5hdGkub3JpZy5ib290LnBlcmYsIE1lYXN1cmUsIFZhbHVlLCBEeHk6YFM6cGAsIGZhY3Rvcl9rZXkgPSBUUlVFKQ0KcDwtZ2dwbG90KGJyaW5hdGkub3JpZy5ib290LnBlcmYgJT4lIGZpbHRlcihNZWFzdXJlICVpbiUgbWVhc3VyZXMpLCBhZXMoeD1NZWFzdXJlLCB5PVZhbHVlLCAgY29sb3I9TWVhc3VyZSkpICsgZ2VvbV9ib3hwbG90KCkgKw0KICBnZW9tX3BvaW50KGRhdGE9ZGF0YS5mcmFtZShNZWFzdXJlID0gbWVhc3VyZXMsIFZhbHVlID0gYygxLCAwLCAwLCAxKSksIHNpemU9MywgY29sb3I9J2JsdWUnKSArIHRoZW1lX2J3KCkNCnA8LUFwcGx5RmlndXJlVGhlbWUocCkNClNhdmVTdGRTcXVhcmVGaWd1cmUocCwgJ2ZpZ3Mvb3JpZy1tb2RlbC1ib290c3RyYXAtZXZhbHVhdGlvbi5wbmcnKQ0KcA0KYGBgDQoNCiMgQXNzZXNzIHBlcmZvcm1hbmNlIGluIGV4dGVybmFsIHZhbGlkYXRpb24NCg0KSUkuIExvb2sgYXQgcGVyZm9ybWFuY2UgYXMgYSBmdW5jdGlvbiBvZiBwcmV2YWxlbmNlDQpgYGB7cn0NCnByZXZhbGVuY2VzID0gc2VxKDAuMDEsIDAuOTksIDAuMDEpDQpicmluYXRpLnN5bi5mYWN0b3J5ID0gZ2VuZXJhdGVTeW50aGV0aWNEYXRhRmFjdG9yeShiYXNlREYgPSBicmluYXRpX2NvdmlkX3N0dWR5X3YyLmltcHV0ZWQsIG1ldGhvZCA9ICdib290JykNCmRmc19wcmV2cyA8LSBsYXBwbHkocHJldmFsZW5jZXMsIGZ1bmN0aW9uIChwcmV2KSBicmluYXRpLnN5bi5mYWN0b3J5KHByZXYgPSBwcmV2KSkNCmRmcy5wcmV2cy5wZXJmIDwtIGxhcHBseShkZnNfcHJldnMsIGZ1bmN0aW9uKGRmKSB7DQogIGFzc2Vzc1BlcmYocHJlZGljdGVkID0gcGxvZ2lzKHByZWRpY3QoZ2xtLm9yaWcuZml0LCBuZXdkYXRhPWFzLmRhdGEuZnJhbWUoZGYpKSksIG9ic2VydmVkID0gZGYkU1dBQikNCn0pDQpgYGANCg0KQ2FsY3VsYXRlIHRoZSBkaWZmZXJlbnQgdmFsaWRhdGlvbiBtZWFzdXJlcw0KYGBge3J9DQoNCmRmcy5wcmV2cy5wZXJmLmRmIDwtIGRhdGEuZnJhbWUobWF0cml4KHVubGlzdChkZnMucHJldnMucGVyZiksIG5yb3cgPSBsZW5ndGgocHJldmFsZW5jZXMpLCBieXJvdyA9IFRSVUUpKQ0KY29sbmFtZXMoZGZzLnByZXZzLnBlcmYuZGYpIDwtIG5hbWVzKGRmcy5wcmV2cy5wZXJmW1sxXV0pDQpkZnMucHJldnMucGVyZi5kZiRQcmV2YWxlbmNlID0gcHJldmFsZW5jZXMNCmRmcy5wcmV2cy5wZXJmLmRmIDwtIGdhdGhlcihkZnMucHJldnMucGVyZi5kZiwgTWVhc3VyZSwgVmFsdWUsIER4eTpgUzpwYCwgZmFjdG9yX2tleSA9IFRSVUUpDQpgYGANCg0KUGxvdCB0aGUgZGlzdHJpYnV0aW9uIG9mIHZhbGlkYXRpb24gbWVhc3VyZXMNCmBgYHtyfQ0KbWVhc3VyZXMgPSBjKCdDIChST0MpJywgJ0JyaWVyJywgJ0ludGVyY2VwdCcsICdTbG9wZScpDQpwPC1nZ3Bsb3QoZGZzLnByZXZzLnBlcmYuZGYgJT4lIGZpbHRlcihNZWFzdXJlICVpbiUgbWVhc3VyZXMpLCBhZXMoeD1NZWFzdXJlLCB5PVZhbHVlLCAgY29sb3I9TWVhc3VyZSkpICsgZ2VvbV9ib3hwbG90KCkgKw0KICBnZW9tX3BvaW50KGRhdGE9ZGF0YS5mcmFtZShNZWFzdXJlID0gbWVhc3VyZXMsIFZhbHVlID0gYygxLCAwLCAwLCAxKSksIHNpemU9MywgY29sb3I9J2JsdWUnKSArIHRoZW1lX2J3KCkNCnA8LUFwcGx5RmlndXJlVGhlbWVMYXJnZUZvbnRPbmx5KHApDQpTYXZlTWVkUmVjdEZpZ3VyZShwLCAnZmlncy9jaGFuZ2UtaW4tbWV0cmljcy1kdWUtdG8tcHJldmFsZW5jZS1ib3hwbG90LnBuZycpDQpwDQpwPC1nZ3Bsb3QoZGZzLnByZXZzLnBlcmYuZGYgJT4lIGZpbHRlcihNZWFzdXJlICVpbiUgbWVhc3VyZXMpLCBhZXMoeD1QcmV2YWxlbmNlLCB5PVZhbHVlLCAgY29sb3I9TWVhc3VyZSwgZ3JvdXA9TWVhc3VyZSkpICsgZ2VvbV9saW5lKHNpemU9MS4zKSArIHRoZW1lX2J3KCkNCnA8LUFwcGx5RmlndXJlVGhlbWVMYXJnZUZvbnRPbmx5KHApDQpTYXZlTWVkUmVjdEZpZ3VyZShwLCAnZmlncy9jaGFuZ2UtaW4tbWV0cmljcy1kdWUtdG8tcHJldmFsZW5jZS1saW5lYXJwbG90LnBuZycpDQpwDQpgYGANCg0KIyBDYWxjdWxhdGUgRGVjaXNpb24gQ3VydmUNCmBgYHtyfQ0KY21zLmN1dG9mZnMgPC0gbGFwcGx5KHNlcSgwLjAxLDAuOTksMC4wMSksIGZ1bmN0aW9uKGN1dG9mZikgew0KICByZXQ9Y29uZnVzaW9uTWF0cml4KHByZWRpY3RlZCwgb2JzZXJ2ZWQsIGN1dG9mZj1jdXRvZmYpDQogIHJldCRjdXRvZmYgPSBjdXRvZmYNCiAgcmV0DQogIH0pDQoNCmNhbGN1bGF0ZV9uZXRfYmVuZWZpdF9hdF90aHJlc2hvbGRfZm9yX3NlbnNpdGl2aXR5X2FuZF9zcGVjaWZpY2l0eSA8LSBmdW5jdGlvbihzZW5zLCBzcGVjLCB0aHJlc2hvbGQpIHsNCiAgc2VucyAtICgxLXNwZWMpKih0aHJlc2hvbGQvKDEtdGhyZXNob2xkKSkNCn0NCm5ldGJlbmVmaXRfYWNyb3NzX3RocmVzaG9sZHMgPSBzYXBwbHkoY21zLmN1dG9mZnMsIGZ1bmN0aW9uKGN1dG9mZl92YWx1ZXMpIGNhbGN1bGF0ZV9uZXRfYmVuZWZpdF9hdF90aHJlc2hvbGRfZm9yX3NlbnNpdGl2aXR5X2FuZF9zcGVjaWZpY2l0eShjdXRvZmZfdmFsdWVzJHNlbnMsIGN1dG9mZl92YWx1ZXMkc3BlYywgY3V0b2ZmX3ZhbHVlcyRjdXRvZmYpKQ0KbmV0YmVuZWZpdF92c190aHJlc2hvbGRzID0gZGF0YS5mcmFtZSh0aHJlc2hvbGRzID0gc2FwcGx5KGNtcy5jdXRvZmZzLCBmdW5jdGlvbihjdXRvZmZfdmFsdWVzKSBjdXRvZmZfdmFsdWVzJGN1dG9mZiksIG5ldGJlbmVmaXQgPSBuZXRiZW5lZml0X2Fjcm9zc190aHJlc2hvbGRzKQ0KZ2dwbG90KG5ldGJlbmVmaXRfdnNfdGhyZXNob2xkcywgYWVzKHg9dGhyZXNob2xkcywgeT1uZXRiZW5lZml0KSkgKyBnZW9tX2xpbmUoKSArIHlsYWIoIk5ldCBCZW5lZml0IikgKyB4bGFiKCJUaHJlc2hvbGQgdG8gQWN0IikgKyB0aGVtZV9idygpDQpgYGANCg0K